Highly localized clustering states in a granular gas driven by a vibrating wall 
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An ensemble of inelastically colliding grains driven by a vibrating wall in 2D exhibits density 
clustering. Working in the limit of nearly elastic collisions and employing granular hydrodynamics, 
we predict, by a marginal stability analysis, a spontaneous symmetry breaking of the extended 
clustering state (ECS). 2D steady-state solutions found numerically describe localized clustering 
states (LCSs). A time-dependent granular hydrodynamic simulation shows that LCSs can develop 
from natural initial conditions. The predicted instability should be observable in experiment. 
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Granular flows exhibit fascinating non-equilibrium 
phenomena and continue to attract much interest [jl],^). 
We will concentrate here on the striking tendency of gran- 
ular "gases" (rapid granular flows) to form dense clusters 
H . Clustering results from energy losses by inelastic col- 
lisions, and it is a manifestation of thermal condensation 
instability fjj. Since the discovery of the clustering in- 
stability, the validity of granular hydrodynamics || has 
been under scrutiny. In a freely cooling granular gas, all 
grains eventually come to rest, making the hydrodynamic 
(and even kinetic) description problematic. In a driven 
granular gas hydrodynamics can be conveniently tested 
on its steady states. The simplest system of this kind 
is a submonolayer of grains in 2D, driven by a vibrating 
side wall at zero gravity. This and related "test bed" 
systems were investigated by molecular dynamic (MD) 
simulations and in experiment For sufficiently 

high average densities, ECS, was observed in these works 
apart from the driving wall. The basic physics of the ECS 
is simple. Because of the inelastic collisions the granular 
temperature decreases with the distance from the driving 
wall. To maintain the momentum balance, the granular 
density increases with this distance, reaching the maxi- 
mum at the opposite ("elastic") wall. When the density 
contrast is large enough, the enhanced density region is 
observed as the ECS. 

Comparisons of the steady-state density profiles ob- 
tained in MD simulations of this problem with those pre- 
dicted by granular hydrodynamics showed that hydrody- 
namics is valid only in the limit of nearly elastic collisions 
|B-BLB. This limit has not been fully explored, and it 
is non-trivial. To show it, we perform a stability analy- 
sis of this simple granular system. This analysis reveals 
a symmetry-breaking instability of the ECS, and forma- 
tion of LCSs. This result puts this system in the list of 
pattern- forming systems out of equilibrium jllj . 

Consider a big ensemble of identical spherical grains of 
diameter d and mass m g = 1 rolling on a smooth horizon- 
tal surface of a rectangular box with dimensions L x x L y . 
The number density of grains is n(x, y). For a submono- 
layer coverage the maximum value of n corresponds to the 



(hexagonal) close-packing value n c — 2/{\/2>d 2 ). Three 
of the walls are immobile, and grain collisions with them 
are assumed elastic. The fourth wall (located at x = L x ) 
supplies energy to the granulate. We will consider two 
different models of energy supply, see below. The en- 
ergy is being lost through inelastic hard-core grain col- 
lisions. We neglect the grain rotation and parameterize 
the inelasticity of grain collisions by a constant normal 
restitution coefficient r. 

We assume a strong inequality 1 — r 2 <C 1, which 
makes a hydrodynamic description valid There- 
fore, steady states of the system can be described by the 
equations of momentum and energy balance: 



p = const , V • (k VT) = I . 



(1) 



where p is the granular pressure, k is the thermal con- 
ductivity, / is the rate of energy losses by collisions and 
T is the granular temperature. To proceed, one needs an 
equation of state p = p (n, T) and relations for n and / 
in terms of n and T. In the low-density limit, n <C n c , 
these relations can be derived from the Boltzmann equa- 
tion The high-density limit, n c — n -C n c , was con- 
sidered by Grossman et al. || . They also suggested con- 
venient interpolations between the low- and high-density 
limits, and verified them by a detailed comparison with 
MD simulations. We will adopt this practical approach 
(see, however, Ref. 0]). In our notation 



p = nT ■ 



(2) 



K = (n/l) n (al + d) 2 T 1 / 2 and I = ( M / 7 1) (1 - r 2 ) n T 3 / 2 . 
Here I is the mean free path of the grains, 



/ = 



V8nd n c - an 



(3) 



a = 1 — (3/8) 1 ^ 2 , and a and 7 are numerical factors of 
order unity. Grossman et al. Q found that a ~ 1.15 and 
7 ~ 2.26. The value of /1, another numerical factor of 
order unity, is irrelevant in the steady-state problem. 

The boundary conditions include the no-flux condi- 
tions V„T = at the "elastic" walls x — 0, y = 



1 



and y = L y . Previously, the "thermal" wall condition 



T = const was used at x = L x |6|-|qjl3|. We will use a 
different boundary condition, to simulate the vibrating 
wall more directly. Our main results, however, will be 
shown to hold for the thermal wall as well. 

The problem of computing the energy flux q from 
a vibrating wall into granulate has been addressed in 
several works Q. Let the wall oscillate sinusoidally: 
x = L x + A cos uit. For small area fractions the granulate 
near this wall is in the dilute limit. We assume A -C I, so 
the vibrating wall does not generate any collective mo- 
tions in the granulate. Grain collisions with the vibrating 
wall are assumed elastic. Also, ui is much larger than the 
rate of granular collisions near the vibrating wall, T 1 / 2 /I, 
so there are no correlations between two successive grain 
collisions with the wall. The limit Aoj <C T 1 / 2 was con- 
sidered by Kumaran |jl5|| for a non-zero gravity. A direct 
calculation analogous to his, but for zero gravity, yields 
q = (2/tt) 1 / 2 A 2 ^nT 1 / 2 |6|. In the language of hydro- 
dynamics, q is the heat flux at the wall: 



ndT/dx = q at x = L x 



(4) 



Finally, f Q x J Q v n(x, y) dxdy = {n) L x L y is normaliza- 
tion condition, where (n) is the average grain density. 

Using Eq. (^), we eliminate T in favor of n and p. 
In its turn, p can be eliminated by integrating Eq. (|l|) 
over the whole box and using the Gauss theorem and Eq. 
(|J) . It is convenient to write the governing equations in a 
scaled form. Introduce scaled coordinates: r/L. x — > r so 
that the box dimensions become 1 x A, where A = L y /L x 
is the aspect ratio of the box. Introducing the (scaled) 
inverse granular density z (x, y) — n c /n (x, y), we obtain 



V • [F{z)Vz) = CQ(z). 
The boundary conditions are 

V„z = at x = , y = and y = A . 

and 



(5) 



(6) 



G(z 



dz 
dx 



1 A 

/ / Qdxdy 



(7) 



J H[z(l,y)]dy 
o 



The normalization condition becomes 



o 



z 



-1 



dx dy = f A . 



while functions F, G, H and Q are the following: 

(z 2 + 2z - \)[az(z 



F(z) 



1 



732/3(2 - a)}' 



(z-a)(z-iy/ 2 z 3 / 2 {z + l) 5 / 2 



(8) 



(9) 



G{z) 



{z 2 + 2z- l)[az{z - 1) + v/32/3(z - a)f 
z(z - a)(z - l)(z + l) 2 



(10) 



, s F(z) , s (z-a)(z- l) 1 / 2 . , 

H ( z ) = tttt and Q( z ) = , i %m m ■ ( n ) 



G(z) 



[Z + l)'V2 z l/2 



Finally, C = (32/37) {L x /d) 2 (1-r 2 ). The other two gov- 
erning parameters are the grain area fraction / = (n) /n c 
and A. Notice that the steady-state density distributions 
are independent of A and u>. 

Eqs. (||)-((|) make a closed set. Their ID (y- 
independent) solution z — Z(x) is described by equations 



(FZ')' =CQ, Z' 



x=0 



and / Z' x dx = / , (12) 



where primes stand for the ^-derivatives. Eq. (0) is now 
satisfied automatically. Eqs. (|12j ) coincide with those 
obtained by Grossman et al. |g] for a thermal wall at 
x = 1. Therefore, the density profiles of the ID states 
coincide for the different types of driving. Eqs. jl^ ) 
can be solved analytically in the high- and low-density 
limits H . These solutions clearly show that criterion fj(i[| 
for validity of hydrodynamics is equivalent to a strong 
inequality 1-r 2 « 1. 

Most interesting among the ID states is the state with 
a dense cluster (ECS) located at the elastic wall x = 0, 
and a low-density region elsewhere. In this case Eqs. ( |T^ ) 
should be solved numerically. Examples are presented in 
Ref. H , and a similar clustering state (CS) was observed 
experimentally ||] . The main objective of this Letter is to 
show that this state (uniform in the y-direction) can give 
way, via a spontaneous symmetry breaking, to CSs highly 
localized in the y-direction. First, a marginal stability 
(MS) analysis will show, in some region of parameters, 
loss of stability of the ECS . Then, solving Eqs. (§)-(§) 
numerically, we will find LCSs. Finally, a time-dependent 
granular hydrodynamic simulation will show that highly 
localized CSs can develop from natural initial conditions. 

Linearizing Eqs. (§)-(§) around the ECS z — Z(x) 
and looking for a small correction in the form of 
ip m {x) cos (7rrny/A), where m — 1, 2, . . . , we obtain 



F4>" ~(CQ 2 



+ n 2 m 2 A- 2 F) 



0. 



(13) 



where <fi — Fip m and index Z means the z-derivative 
evaluated at z = Z(x). The boundary conditions are 

0'U=o = O, {FGfi + Z l (FGz-GF z ),]>)\ x ^ 1 =Q. (14) 

Functions F and G which enter Eqs. (|l3|) and (Q) are 
evaluated at z = Z(x). 

For fixed values of £, / and m, Eqs. (O) and ( p^ ) 
represent a linear eigenvalue problem for A , or simply 



2 



for the aspect ratio A. We denote these eigenvalues by 
A m , to = 1, 2 .... Obviously, A m (C, f) = mA^CJ). 
The eigenvalues A m (when they exist) represent the crit- 
ical values of the aspect ratio: at A > A m the ECS 
looses stability. For the mode to = 1 (or A/2, one half 
of the wavelength across the system in the y-direction) 
the critical value Ai is the lowest. Figure 1 shows, for 
different values of £, the MS curves A = Ai(/) that 
we found numerically. Above the MS curves LCSs must 
develop. Interestingly, the ECS remains linearly stable 
for any A beyond a finite interval of the area fractions 
/i(£) < f < f 2 (C) such that /i > and f 2 < 1. As £ 
increases, the instability interval (/i,/2) shrinks, while 



the minimum value of Ai decreases: A^' 



(mm) 



uoc- 1 / 2 



|L7| . For sufficiently large C, A J""™- 1 becomes less than 1. 
See Fig. 2 which shows the instability tongues m = 1 and 
to = 2 for C = 5 • 10 4 . For higher to we obtain modes 
3A/2, 2A, 5A/2, . . . which fit in boxes with increasingly 
larger aspect ratios, A > toAi. 
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FIG. 1. Marginal stability curves Ai(/) for different values 
of C. The values of Ai are multiplied by C 1 ^ 2 . 
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FIG. 2. Marginal stability (MS) curves A m (/) for m = 1 
(solid line) and 2 (dashed line) and £ = 5-10 4 . Density profiles 
corresponding to points 1,2,3 and 4 are shown in Fig. 3. The 
dotted line shows the MS curve m — 1 for the "thermal wall" . 

When / <C min (1, £ -1 / 2 ) , Ai(/) can be found ana- 
lytically. In this case the whole system is in the dilute 
limit, z> 1 (still, it is necessary to account for the sub- 
leading terms). In addition, Z(l) - Z(0) < Z(0) in this 
case, so Taylor expansion of Z(x) and ip( x ) U P to x suf- 



fice. After some algebra, Eqs. (§J) - @ yield 

/£ 2 / 4 (l + a)£f 3 Y 1/2 
Al=7r (l^" J (15) 

It follows from Eq. that /i(£) = 3q 2 (1 + a) £-\ 

The same instability appears when the wall x = 1 (in 
scaled units) is "thermal" . Solving the corresponding 
eigenvalue problem [where the second boundary condi- 
tion in Eq. ( |l4| ) is replaced by <p(x = 1) = 0], we obtained 
instability tongues similar to those for the vibrating wall, 
but more narrow. Figure 2 shows the instability tongue 
to = 1 for C = 5 • 10 4 . Noticeable is the coincidence of 
the to = 1 curves at intermediate / for the two types 
of driving. It results from a strong localization of the 
eigenfunction <f>(x) near the elastic wall x = at large 
£ and intermediate /. The exact form of the boundary 
condition at x = 1 is not important in this regime. Fi- 
nally, for the thermal wall the LCS is stable for any A if 
/ <C min (1, £ -1 / 2 ), in contrast to the vibrating wall. 

In the rest of the paper we will deal with the vibrat- 
ing wall. Within the instability tongues of Fig. 1 the 
MS analysis is invalid. Besides, it can miss a subcritical 
bifurcation outside of the instability tongues. Therefore, 
we solved the 2D steady-state equations (||)-(||) numer- 
ically (a nonlinear Poisson solver, Newton's iterations), 
exploring some parts of the parameter plane (/, A) of 
Fig. 2. Figure 3 shows the density profiles of 4 typical 
steady states with to = 1 and 2. Highly localized nonlin- 
ear A/2- and A-states are evident in Fig. 3a and b. The 
maximum/minimum density ratio along the elastic wall 
x = reaches about 21 in these examples. 



FIG. 3. Steady-state density profiles (gray scale, separate 
for each picture) corresponding to points 1 (a), 2 (b), 3 (c) 
and 4 (d) of Fig. 2. The maximum (minimum) density at 
the wall x = is 0.76 (0.036) (a and b), 0.48 (0.21) (c) and 
0.54 (0.10) (d). The gas density at the vibrating wall x = 1 
is close to 4 • 10 -3 for all profiles. 

In general, we found that when crossing the MS curve 
A = Ai(f) from the left (along the line A = 2), or from 
below, one goes continuously from an ECS to a "weakly 
2D" A/2-state. This implies a supercritical bifurcation. 
However, when moving from the right to the left along 
the line A = 2, nonlinear A/2- and A-states appear inside 
the linear stability regions of the ECS and of the mode 
to = 1, respectively, and coexist with the ECS (with the 
mode to = 1, respectively). These findings give evidence 
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for bistability and subcritical bifurcations. Examples of 
subcritical A/2- and A-states are shown in Fig. 3c and 
d. Super- and subcritical LCSs were also observed for 
A = 3. 

A mirror reflection of Fig. 3a with respect to y = 
makes A = 2 and produces a A-state similar to Fig. 2b. 
Furthermore, extending Fig. 2b periodically in the y- 
direction, we obtain a periodic chain of LCSs which fit in 
boxes with increasingly larger aspect ratios. When the 
aspect ratio goes to infinity, the periodic cluster chain 
becomes infinite. Cluster chains with different periods 
can fit in boxes with large enough aspect ratios, therefore, 
an interesting selection problem appears, like in other 
pattern-forming systems . 

We performed a series of time- dependent hydrody- 
namic simulations (with A = 1,2 and 3) which showed 
that the LCSs are dynamically stable and develop from 
natural initial conditions. We will briefly report here a 
single simulation with A = 2. The full hydrodynamic 
equations were solved with the same constitutive rela- 
tions and boundary conditions as in the steady state 
analysis. Instead of the shear viscosity in the momen- 
tum equation we accounted for a small rolling friction 
force — nv/r. An extended version of the compressible 
hydro code VULCAN @ was used. 

The initial scaled density in this example was 
n(x,y,t = 0) = / + 0.1/ cos(7ry) (independent of x). 
Figure 4 shows the density evolution. A cluster develops 
near the elastic wall x = 0. With time it becomes local- 
ized in the y-direction and approaches the steady-state 
profile shown in Fig. 3b. 
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FIG. 4. Density evolution for £ = 5 • 10 4 , A = 2 and 
/ = 0.0235. Shown are the density profiles (gray scale, sepa- 
rate for each picture) at scaled times 100 (a), 500 (b), 1,000 
(c) and 1, 290(d). The maximum (minimum) density at the 
wall x = is 0.25 (0.14) (a), 0.46 (0.072) (b), 0.66 (0.040) (c) 
and 0.74 (0.036) (d). The gas density at the vibrating wall 
x = 1 is close to 4 ■ 10 -3 for all profiles. 

In summary, we predict a spontaneous transition from 
an extended to highly localized clustering states in a 
driven submonolayer granular system. The transition 
should occur when the aspect ratio is large enough (see 
Fig. 1). It is insensitive to the vibration frequency and 
amplitude, and only weakly depends on the type of the 
driving wall. To observe this transition in experiment, 
one should minimize the role of the rolling/sliding fric- 



tion M, unaccounted for in our model. The frictional 
energy losses are proportional to T, while the collisional 
energy losses are proportional to T 3 / 2 . Therefore, one 
should work with higher granular temperatures (that is, 
larger Aid). 

When 1 — r 2 is not small, the normal stress difference, 
non-Gaussianity in the velocity distribution and possible 
lack of scale separation become important. The role of 
these effects in the symmetry-breaking instability should 
be the subject of further studies. 

Finally, the aspect ratios used in Refs. were al- 

ways lower than the threshold values for the instability 
A[ mm '. As the result, the instability was suppressed by 
granular heat conduction in the y direction. 
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